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The highly- frustrated spin-| quantum Heisenberg model with both nearest (Ji) and next-nearest 
(J2) neighbor exchange interactions is revisited by using an extended variational space of projected 
wave functions that are optimized with state-of-the-art methods. Competition between modulated 
valence-bond crystals (VBCs) proposed in the literature and the Dirac spin liquid (DSL) is investi- 
gated. We find that the addition of a small ferromagnetic next- nearest-neighbor exchange coupling 
I J2I > 0.09Ji leads to stabilization of a 36-site unit cell VBC, although the DSL remains a local min- 
imum of the variational parameter landscape. This implies that the VBC is not trivially connected 
to the DSL; instead it possesses a non-trivial flux pattern and large dimerization. 

PACS numbers: 75.10.Kt, 75.10.Jm, 75.40.Mg 



Introduction. It is well known that, on the nonbipar- 
tite two-dimensional Kagome lattice, the combination of 
low spin {S = 1/2), low coordination number {z = 4), 
and frustrating anti-ferromagnetic (AF) exchange inter- 
actions lead to extremely strong quantum fluctuations. 
It is, however, a widely debated and long-standing the- 
oretical issue whether the ground state of the nearest- 
neighbor (n.n.) spin-1/2 quantum Heisenberg antiferro- 
magnet (QHAF) on the Kagome lattice is a spin dis- 
ordered state (quantum spin liquid), ^'^ which preserves 
spin rotation and lattice space group symmetry, or in- 
stead a valence bond crystal (VBC),^ which breaks 
lattice symmetries. On the experimental side, studies 
on a nearly perfect spin-1/2 Kagome compound, Her- 
bertsmithite ZnCu3(OII)6Cl2,^~^^ reveal the absence of 
any spin ordering down to 50 mK despite a sizable n.n. 
AF exchange coupling (J ^ 180K) between spin-1/2 
moments of Cu^^. In particular, Raman spectroscopic 
data on Herbert smithite points towards a gapless, alge- 
braic spin liquid state. This lends support to the view 
that the ground state of the n.n. spin-1/2 QHAF model 
on the Kagome lattice is a long-range resonating- valence 
bond state. Within a class of variational projected wave 
functions, a particular gapless spin liquid belonging to 
the class of algebraic spin liquids, the U(l) Dirac state, 
has been claimed to possess the lowest energy. ^^'^^ In 
such a state, the (mean-field) Fermi surface collapses 
to two points, where the spectrum becomes relativistic 
with Dirac conical excitations. On the contrary, a recent 
study of the n.n. spin-1/2 QHAF model using density- 
matrix renormalization group (DMRG),^^ establishes the 
ground state to be a singlet-gapped spin liquid, suppos- 
edly with a Z2 low-energy gauge structure. 

From the experimental point of view, the weak fer- 
romagnetism observed in Herbertsmithite has been at- 
tributed to the ferromagnetic (FM) nature of the next- 
nearest-neighbor (n.n.n.) coupling between Cu^^ ions 
in the Kagome layers. This model was investigated in 



Ref. [20] by using projected wave functions and it has 
been found that, above a certain critical n.n.n. FM cou- 
pling, a gapless spin liquid with a large circular spinon 
Fermi surface, named here as uniform projected Fermi 
sea (PFS), is stabilized as opposed to the U{1) Dirac 
state. Furthermore, this state undergoes a small dimer- 
ization, which lowers slightly its energy. The same n.n.n. 
FM model was also recently investigated using a quantum 
dimer model approach in Ref. [21], showing consistently 
that a 36-site cell VBC order is favored. 

In this Rapid Communication, we revisit the spin-1/2 
QHAF with the inclusion of n.n.n. exchange interac- 
tions using a more extended variational space of pro- 
jected wave functions that may be optimized by using 
the technique of Ref. [22]. In the following, we will limit 
to non-magnetic variational states and, therefore, we will 
not consider possible instabilities toward magnetically or- 
dered states. Our main result is that the addition of 
ferromagnetic n.n.n. exchange coupling leads to the sta- 
bilization of a 36-site unit cell VBC (in agreement with 
the results of Ref. [21]) over an extended ferromagnetic 
region which starts from a very weak coupling. Although 
being a dimerization of the uniform PFS, our VBC does 
not arise from a local instability of the latter. In other 
words, it is not trivially connected to it in the varia- 
tional parameter landscape; instead it possesses a non- 
trivial flux pattern and dimerization. Moreover, we find 
that the level crossing between the PFS and the U{1) 
Dirac state (once suitably extended with n.n.n. hopping) 
occurs at nearly twice the value previously reported. 
For AF n.n.n. exchange coupling, the inclusion of n.n.n. 
hoppings in the U{1) Dirac state leads to a considerable 
lowering in energy, which becomes more pronounced with 
increasing the AF coupling. Moreover, no VBC order is 
found in the AF n.n.n. coupling region. 

Model and wave function. The Hamiltonian for spin- 
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1/2 quantum Heisenberg J1—J2 model is 



H 



(ij) 



J2 ^ Si 

{(ij)) 



(1) 



where {ij) and ((^j)) denote sums over n.n. and n.n.n. 
neighbor sites, respectively. In the fohowing, we wih con- 
sider Ji > and both FM and AF superexchange J2; all 
energies will be given in units of Ji. 

The variational wave functions are defined by project- 
ing noncorrelated fermionic states: 



(2) 



where Vq = Hill ~ ^i,t^i,i) ^^^^ Gutzwiller 

projector enforcing the one fermion per site constraint. 
Here, |^MF(Xij7 ^ij? m)) is the ground state of mean-field 
Hamiltonian containing chemical potential, hopping, and 
pairing terms: 



In this work, all states that we consider include hop- 
ping terms only, i.e., Xij up to second neighbors. The 
effect of including BCS pairing terms is discussed at the 
end of this Rapid Communication. Cases in which the 
translational symmetry is explicitly broken will also be 
considered, so as to include VBC states. 

Different spin-liquid and VBC phases correspond to 
different patterns of distribution of Xij ^^^d A^j on the 
lattice bonds; they are the ansatz of a given state and 
serve as the variational parameters that are optimized 
within the variational Monte Carlo scheme to find the 
energetically best state. It is worth mentioning that 
this method allows us to obtain an extremely accurate 
determination of variational parameters. All parameters 
belonging to one class (i.e., with the same magnitude) 
are generically labelled as xa- 

Results. We have performed our variational calcula- 
tions on a 576-site (i.e., 36 x 4 x 4) cluster with mixed 
periodic-antiperiodic boundary conditions. Such a clus- 
ter accommodates all possible VBC supercells proposed 
in the literature. In addition, it ensures non-degenerate 
wave functions at half-filling. 

For n.n. spin- 1/2 QHAF, among the class of 
n.n. translationally symmetric, non-chiral, gapless spin 
liquids, the U{1) Dirac state is given by the ansatz in 
Fig. 1(a). Due to fiux cp being and tt (exp (z(/p) = 
ripiaquetteXA) through the triangles and hexagons re- 
spectively, it is denoted as [0,7r]. Its energy per site is 
E/Ji = -0.42866(1). The n.n. uniform PFS state has 
no fiux through any plaquette and is therefore denoted 
as [0,0], its energy per site is E/Ji = -0.41197(1).^^'^^ 
In this work we study only gapless states in particular 
those with a /7(1) low energy gauge structure. However, 
we believe that by performing a case by case projected 
wave function study of all possible (a few hundred) Z2 
spin liquids on the Kagome lattice, one can identify vari- 
ationally the state found in Ref. [19] using DMRG. 
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FIG. 1. The U{1) DSL ansatz, solid (dashed) bonds denote 
positive (negative) hoppings (a). The unit cell is doubled to 
accommodate [0, tt] flux. Cases with n.n.n. hopping are also 
reported: the only possible (non-chiral) space group symmet- 
ric states built upon [0, tt] state are [0, tt; tt, 0] (b) or [0, tt; 0, tt] 
(d) and those upon the uniform [0,0] state (i.e., PFS) are 
[0,0;7r,7r] (c) or [0,0; 0,0] (e). 



With the aim of investigating the effect of an additional 
n.n.n. exchange coupling (of both AF and FM type), 
we first extend the [0, tt] and [0, 0] states. While previ- 
ous studies, considered wave functions with n.n. cou- 
plings only, here we include in addition n.n.n. couplings 
in the mean-field Hamiltonian (3), which also leads to 
space group symmetric, non-chiral, gapless spin liquids, 
see Fig. 1. Two new plaquettes (triangles abc and acd 
in Fig. 1) appear upon the inclusion of n.n.n. couplings. 
Space group symmetric, non-chiral spin liquids may now 
be labelled by four fiuxes (but only three are indepen- 
dent) (i.e., by [a, /3; 7, ^]: a and f3 are fiuxes through orig- 
inal triangles and hexagons, respectively; 7 and 5 instead 
are fiuxes through triangles abc and acd^ respectively). 
The only possible states built upon the [0,7r] state are 
[0, tt; tt, 0] or [0, tt; 0, tt] and those upon the [0, 0] state are 
[0,0;7r,7r] or [0,0; 0,0] (see Fig. 1). Notice that for both 
DSL and PFS, the two states with different 7 and S fiuxes 
are related by a change of sign in Xn.n.n.- The energeti- 
cally lower states will depend upon the actual value of the 
ratio J2/ Ji- This extension does not modify the topolog- 
ical properties associated with the wave functions, such 
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FIG. 2. (a) Optimized Xn.n.n./Xn.n. vs J2 for the extended 
PFS and DSL states of Fig. 1. A zoom around J2 = is 
shown in (b). The optimized n.n. (c) and n.n.n. (d) hopping 
parameters vs J2 for the 36-site supercell VBC are also re- 
ported, xi = 1 is the reference bond. Only xii < 0, which 
implies [0, 0; tt, tt] flux in the P hexagons of Fig. 4. 



as the Dirac spectrum and the large spinon Fermi sur- 
face. Most importantly, the inclusion of n.n.n. hopping 
parameters leads to lowering of the variational energies, 
via an optimal tuning of Xn.n.n. as a function of J2 [see 
Fig. 2(a)]. It is important to note that we purposely re- 
strict our calculations to small enough J2/J1, since for 
larger n.n.n. couplings (of both AF and FM type), it is 
probable that Neel states are energetically favored, and 
consequently our treatment becomes insufficient. 

Point D in Figs. 2(b) and 3 marks a transition between 
the [0, 0; 0, 0] and [0, 0; tt, tt] states, and point E, the tran- 
sition between [0, tt; 0, tt] and [0, tt; tt, 0] states, both oc- 
curring at J 2 7^ 0. Therefore, we find a finite Xn.n.n. 
even for the n.n. spin- 1/2 QHAF [see points F and G in 
Fig. 2(b)]. We mention that these extended wave func- 
tions with n.n.n. hoppings lead to slightly lower energies, 
namely E / Ji = -0.42872(1) for the [0,7r;0,7r] state and 
E/Ji = -0.41209(1) for the [0,0;7r,7r] state. 

Due to negative n.n.n. spin-spin correlations of [0, tt] 
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FIG. 3. Energy vs J2 for spin liquid (see Fig. 1) and VBC 
states (see Fig. 4(c) and (d)). 



state and positive for the [0, 0] state, a level crossing oc- 
curs at J2/J1 ~ —0.16^^ (see point B in Fig. 3). How- 
ever, the addition of the n.n.n. hopping shifts the level 
crossing between the reference spin liquids, the [0, 0; 0, 0] 
and [0,7r;0,7r] states to J2/J1 ^ —0.335 (see point A 
in Fig. 3). The question of global and local instabil- 
ity of these spin-liquid states towards a VBC ordering is 
now thoroughly addressed. In contrast to previous stud- 
17,20 aimed at checking only the local instabili- 

ties of spin liquids towards various dimerization patterns 
via the imposition of a small bond amplitude modula- 
tion (5%-10%) of xa, we make a complete optimization 
of the parameters to detect a possible stabilization of 
VBC states. 

In the 12-site supercell, all bonds connected by Dq op- 
erations have same magnitude [see Fig. 4(a)], leading to 
three classes of different bonds. In the 18-site supercell 
[see Fig. 4(b)], there are only two classes of bonds. In 
the 36-site supercell with n.n. couplings [see Fig. 4(c)], 
there are six classes, given the Dq symmetry about the 
hexagon C. By adding n.n.n. bonds and preserving this 
symmetry, we obtain six more independent bonds [see 
Fig. 4(d)] 

In our analysis, we start from a large number of arbi- 
trary different points (amplitude modulation of xx) the 
variational space and thoroughly scan the landscape. As 
a consequence we find that for the n.n. spin-1/2 QHAF, 
the [0,7r], and [0,0] states are locally and globally sta- 
ble with respect to 12-site,^ IS-site,^'^^ and 36-site^'^ 
supercell dimerizations. Indeed, the optimization pro- 
cedure always gives back the uniform |xa| = 1 state. 
We bring attention to the fact, that the 36-site super- 
cell considered by us has a much larger variational space 
(six different hoppings) as compared to the ones studied 
in literature, which considered a dimerization only along 
the hexagon P in Fig. 4(c). 

Upon the inclusion of a finite FM J2, we detect the 
appearance of another competing state with broken sym- 
metries which is stabilized and is the lowest in energy for 



4 




FIG. 4. 12-site supercell (a), with three different parameters 
for the hopping; 18-site supercell, with two different parame- 
ters; 36-site supercell, with six n.n. parameters (c), and with 
six n.n.n. parameters (d). 



J2 < —0.09 (see point C in Fig. 3). This state is found 
to be a 36-site supercell VBC, shown in Fig. 4(c) and 



4(d). It breaks translational symmetry in the magnitude 
of the n.n. and n.n.n. order parameters and the [7, 6] 
fluxes but preserves the rotation and reflection symme- 
try in both magnitude of xa and [a, /3; 7, S] fluxes. The 
corresponding state has 12 different hopping parameters. 
Although we obtain this VBC as a dimerization of the 
[0, 0; 0, 0] state, it is not a local instability of it. Instead, 
it possesses a large bond amplitude modulation in the 
extended variational space of n.n. and n.n.n. order pa- 
rameters and selects a flux pattern with [7, 6] fluxes be- 
ing [tt, tt] in the P hexagons which form a honeycomb 
lattice, see Fig 4(d). On the contrary, all other hexagons 
have [0, 0] fluxes. The optimized xa as function of J2 are 
shown in Fig. 2(c) and 2(d). In a previous work, which 
investigated the effect of a J2 FM exchange coupling us- 
ing projected variational wave functions, it was found 
that a dimer modulation leads to an energy minimum at 
J2/J1 ^ —0.16, for approximately 4% bond amplitude 
modulation. In contrast, we find a different VBC wave 
function, which is stabilized starting from a very weak 
FM n.n.n. coupling. As mentioned above, this state 
possesses a very large 36-site modulation, leading to a 
significant large gain in energy. 

We finish by considering the case of an AF J2 coupling. 
Our study reveals the absence of symmetry breaking, in- 
stead we find a gapless state with Dirac fermions, the 
[0, tt; tt, 0] state. Upon optimizing Xn.n.n./Xn.n., this state 
gets a significantly lower energy than the n.n. [0, tt] state, 
this gain becoming more pronounced for larger J2 (see 
Fig. 3). Finally, the addition of a BCS pairing term of 
the 5- wave type in the [0, tt; tt, 0] wave function for J2 AF 
is also studied, and our calculations show that such an 
inclusion always increases the energy. However, effect of 
including other forms of pairing terms which might sta- 
bilize a gapped spin liquid or a VBC in the J2 AF model 
is left as direction of future research. This might provide 
a reconciliation with the exact diagonalization results of 
Ref. [24] , which point towards an opening of a gap upon 
addition of a small AF J2 coupling. 

In summary, we investigated the spin- 1/2 QHAF on 
the Kagome lattice by using improved variational wave 
functions. We found that a VBC is stabilized when a 
n.n.n. ferromagnetic superexchange coupling is consid- 
ered. This state possesses a non-trivial distribution of 
hopping parameters and flux pattern. 
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